clear all
set more off
set type double, perm
capture log close

global d = 0

if $d == 0 { 
global topdir "YOUR DIRECTORY"
global datadir "${topdir}/data"
global outdir "${topdir}/output"
cd "${outdir}"
}

use group weight aafqt using "${datadir}/nlsy_final.dta", clear

// log 
log using bstest_afqt_distribution.log, replace

// set seed
set seed 20220331

// record original sample size 
count
global n = r(N)

// expand by sample weight
gen weight_simple = int(weight/100)
expand weight_simple


cap program drop bstest
program bstest, rclass // test distribution statistics
	
	** preserve
	preserve

	** distribution statistics
	foreach g in WM79 WM97 WW79 WW97 {
	sum aafqt if group=="`g'", d
	foreach stat in mean sd skewness kurtosis p1 p5 p10 p25 p50 p75 p90 p95 p99 {
	local `stat'_`g' = r(`stat')
	return scalar `stat'_`g' = ``stat'_`g''
	}
	}
	
	** return - difference between nlsy79 and nlsy97
	foreach g in WM WW {
	foreach stat in mean sd skewness kurtosis p1 p5 p10 p25 p50 p75 p90 p95 p99 {
	return scalar d_`stat'_`g' = ``stat'_`g'97' - ``stat'_`g'79'
	}
	}

	** restore
	restore
	
end


bootstrap d_mean_WM=r(d_mean_WM) d_sd_WM=r(d_sd_WM) d_skewness_WM=r(d_skewness_WM) d_kurtosis_WM=r(d_kurtosis_WM) ///
	d_p1_WM=r(d_p1_WM) d_p5_WM=r(d_p5_WM) d_p10_WM=r(d_p10_WM) d_p25_WM=r(d_p25_WM) d_p50_WM=r(d_p50_WM) ///
	d_p75_WM=r(d_p75_WM) d_p90_WM=r(d_p90_WM) d_p95_WM=r(d_p95_WM) d_p99_WM=r(d_p99_WM) ///
	d_mean_WW=r(d_mean_WW) d_sd_WW=r(d_sd_WW) d_skewness_WW=r(d_skewness_WW) d_kurtosis_WW=r(d_kurtosis_WW) ///
	d_p1_WW=r(d_p1_WW) d_p5_WW=r(d_p5_WW) d_p10_WW=r(d_p10_WW) d_p25_WW=r(d_p25_WW) d_p50_WW=r(d_p50_WW) ///
	d_p75_WW=r(d_p75_WW) d_p90_WW=r(d_p90_WW) d_p95_WW=r(d_p95_WW) d_p99_WW=r(d_p99_WW), ///
	nodots size($n) reps(2000) saving(bstest, replace): bstest

	
log close
exit
